#ifndef INITIALIZE
#define INITIALIZE
/// Includes cuda
#include <cutil_inline.h>
/// Includes projects
#include "..\common\svmCommon.h"
#include "svmKernels.h"
#include "svmTrain.h"

__global__ void initF(int nPoints,
					  float* F,
					  float* label)
{
	unsigned int index = threadIdx.x;
	unsigned int stride = blockDim.x;
	while(index < nPoints) {
		F[index] = -label[index];
		index += stride;
	}
}

template<KernelType type> 
__global__ void preComputeKernel(//float* devData, int devDataPitchInFloats,
								 int nPoints,
								 //int nDimension,
								 float parameterA, float parameterB, float parameterC,
								 float* devKernelDiag,
								 float* devSelfDot)
{
	unsigned int index = __mul24(blockIdx.x, blockDim.x) + threadIdx.x;

	if (index < nPoints) {
		if(type == GAUSSIAN) {
			Gaussian kf;
			devKernelDiag[index] = 1.0f;
			devSelfDot[index] = kf.selfDot(index);
		}
		if(type == LINEAR) {
			Linear kf;
			devKernelDiag[index] = kf.selfK(index, parameterA, parameterB, parameterC);
			//devKernelDiag[index] = kf.selfKernel(devData + index, devDataPitchInFloats, devData + (nDimension * devDataPitchInFloats), parameterA, parameterB, parameterC);
		}
		if(type == POLYNOMIAL) {
			Polynomial kf;
			devKernelDiag[index] = kf.selfK(index, parameterA, parameterB, parameterC);
			//devKernelDiag[index] = kf.selfKernel(devData + index, devDataPitchInFloats, devData + (nDimension * devDataPitchInFloats), parameterA, parameterB, parameterC);
		}
		if(type == SIGMOID) {
			Sigmoid kf;
			devKernelDiag[index] = kf.selfK(index, parameterA, parameterB, parameterC);
			//devKernelDiag[index] = kf.selfKernel(devData + index, devDataPitchInFloats, devData + (nDimension * devDataPitchInFloats), parameterA, parameterB, parameterC);
		}
	}
}


/// For the first step, we know alpha1Old == alpha2Old == 0, and we know y1 != y2
/// where labels[iLow] = -1 and labels[iHigh] = 1 -> eta > 0
/// So we just boil down the algebra, any value between zero and C is acceptable
template<class Kernel>
__global__ void takeFirstStep(float* devResult,
							  float* devKernelDiag,
							  //float* devData, int devDataPitchInFloats,
							  //float* devTransposedData, int devTransposedDataPitchInFloats,
							  float* iLowP,
							  float* devAlpha,
							  float cost,
							  int iLow, int iHigh,
							  float parameterA, float parameterB, float parameterC,
							  float* selfDot,
							  Kernel kf)
{ 
	float eta = devKernelDiag[iHigh] + devKernelDiag[iLow];
	//float* pointerB = devData + iLow;
	//float* pointerB = devTransposedData + iLow*devTransposedDataPitchInFloats;
	//float* pointerA = devData + iHigh;
	//float* pointerAEnd = devData + nDimension*devDataPitchInFloats;
	//float phiAB = kf.kernel(pointerA, devDataPitchInFloats, pointerAEnd, pointerB, devDataPitchInFloats, parameterA, parameterB, parameterC);
	//float phiAB = Kernel::K(iHigh, iLow, parameterA, parameterB, parameterC);
	//float phiAB = kf.K(iHigh, pointerB, devDataPitchInFloats, parameterA, parameterB, parameterC, selfDot, iLow);
	//float phiAB = kf.K(iHigh, pointerB, parameterA, parameterB, parameterC, selfDot, iLow);
	float phiAB = kf.K(iHigh, iLowP, parameterA, parameterB, parameterC, selfDot, iLow);
	eta = eta - 2.0f*phiAB;


	float alpha2New = 2.0f/eta; 
	if (alpha2New > cost) {
		alpha2New = cost;
	}

	///alpha1New == alpha2New for the first step
	devAlpha[iLow] = alpha2New;
	devAlpha[iHigh] = alpha2New;

	devResult[0] = -alpha2New; /// label_Low(alpha2_new - alpha2_old)
	devResult[1] = alpha2New; /// label_High(alpha1_new - alpha1_old)
}

#endif
